'''
Parser for the HPeak output files

Description of the file format can be found here:

http://www.sph.umich.edu/csg/qin/HPeak/Readme.html

.hpeak.out:
This is the new main output file. Compared to .allregions.txt, the only difference is an inserted extra column (column 5) 
which indicates the location within the peak that has the highest hypothetical DNA fragment (HDF) coverage). 
Note that the first column only contains the chromosome number, and X and Y are replaced by 23 and 24 for easier 
numerical manipulations.

Fields ---
chrom#: number of chromosome (1 to 24)
chromStart: start point of region
chromEnd: end point of region
length: size of region (equal to chromEnd - chromStart) 

Created on May 20, 2009

@author: mkiyer
'''

def parse_hpeak_file(fhd):
    for thisline in fhd:
        yield parse_hpeak_line(thisline) 
        
def parse_hpeak_line(line):
    '''
    parse a single line and return fields
    '''
    fields = line.strip().split('\t')
    # TODO: what to do with rest of fields?
    return ('chr' + fields[0], int(fields[1]), int(fields[2]))
        
class HPeakAnnotation(object):
    
    def __init__(self):
        pass

    @staticmethod
    def parseFile(fhd):
        for thisline in fhd:
            # ensure validity of this line
            assert thisline
            thisline = thisline.strip()
            assert thisline
            # parse the line
            hpeakline = HPeakAnnotation()
            hpeakline.parseLine(thisline)
            yield hpeakline
    
    def parseLine(self, thisline):
        fields = thisline.strip().split('\t')
        
        chromNum = int(fields[0])
        if chromNum == 23:
            self.chrom = 'chrX'
        elif chromNum == 24:
            self.chrom = 'chrY'
        else:
            self.chrom = 'chr' + fields[0]
        
        self.chromStart = int(fields[1]) - 1
        self.chromEnd = int(fields[2])
        self.peakSummitPos = float(fields[4])
        self.peakHeight = float(fields[5])
        self.negLogPValue = float(fields[6])


class HPeakGeneAnnotation(object):
    
    def __init__(self):
        pass

    def __repr__(self):
        return '<HPeakAnnotation> %s' % str(self.fields)
            
    @staticmethod
    def parseFile(fhd):
        for thisline in fhd:
            if thisline.rstrip().startswith('Chromosome'):
                continue
            a = HPeakAnnotation()
            a.parseLine(thisline)
            yield a
    
    def parseLine(self, line):
        '''
        .annotation.txt:
        
        Fields:
        chrom: chromosome:start-stop (e.g. chr1:50-100, chr2:134897-239487)
        coverage: coverage
        peakHeight: peak height
        TODO: UNKNOWN FIELD: UNKNOWN
        percentGC: GC content of peak
        percentMaskedBases: percent of masked bases
        conservationScore: conservation score
        locationType: location type in genome, one of 'Intergenic', 'Intron', 'Exon', 'TSS region'
        gene1_name: conventional gene name for gene 1
        gene1_genbank_name: genbank acc id of gene 1
        gene1_strand: chromosomal strand of gene1
        gene1_distance: distance from peak to gene 1
        gene2_name: conventional gene name for gene 2, or blank if no 2nd gene exists for this peak
        gene2_genbank_name: genbank acc id of gene 2, or blank if no 2nd gene exists for this peak
        gene2_strand: chromosomal strand of gene 2, or blank if no 2nd gene exists for this peak
        gene2_distance: distance from peak to gene 2, or blank if no 2nd gene exists for this peak        
        '''            
        fields = line.strip().split('\t')
        self.fields = fields
        self.chrom, chrom_interval = fields[0].split(':')
        self.chromStart, self.chromEnd = map(int, chrom_interval.split('-'))
        self.coverage = fields[1]
        self.peakHeight = fields[2]
        # TODO: what is this field?
        self.unknownField3 = fields[3]
        self.percentGC = fields[4]
        self.percentMaskedBases = fields[5]
        self.conservationScore = fields[6]
        self.locationType = fields[7]

        gene_name, gene_genbank_name, gene_strand = fields[8:11]
        gene_distance = int(fields[11])
        self.geneInfo = []
        self.geneInfo.append((gene_name, gene_genbank_name, gene_strand, gene_distance))
        if len(fields[12:16]) == 4:
            gene_name, gene_genbank_name, gene_strand = fields[12:15]
            gene_distance = int(fields[15])
            self.geneInfo.append((gene_name, gene_genbank_name, gene_strand, gene_distance))